### CODE FOR:
### 
### Perceiving Welfare State Sustainability: 
### Fiscal Costs, Group Deservingness, or Ideology?
### 
### VERSION: 16-APR-2024


rm(list=ls())
gc()




# WORKING DIRECTORY -------------------------------------------------------
setwd("C:/Users/miros/Desktop/Research/Deservingness frames")





# LIBRARIES ---------------------------------------------------------------
library(haven)
library(labelled)
library(ggplot2)
library(ggpubr)
library(data.table)
library(sjPlot)
library(grid)
library(gridExtra)
library(stargazer)
library(RColorBrewer)
library(interplot)
library(egg)





# DATA PREPARATIONS -------------------------------------------------------
# __ Importing dataset ----------------------------------------------------
# Note: Data are available for research purposes, including replication, 
#       on the webpage of Norwegian Centre for Research Data (SIKT) at 
#       https://doi.org/10.18712/NSD-NSD2458-V6
SuppA <- read_sav("02-data/NSD2458.sav")

# Removing confusing labels
SuppA <- labelled::remove_labels(SuppA)

# Minor format adjustments
SuppA$id <- as.character(SuppA$id)
SuppA$Periode <- factor(SuppA$Periode, levels = c(1, 2, 3))

# Subsetting only observations from the second wave
# (where experimental design was included)
SuppA <- subset(SuppA, Periode == 2)





# __ Selecting relevant variables -----------------------------------------
SuppA <- SuppA[c(1:2,        #Administrative
                 262:271,    #Experiment
                  65,        #Left-Right
                   7,        #Gender
                   8,        #Age
                   9,        #Education
                  12,        #Personal Gross Income
                  32         #Weight
                 )
               ]





# RESPONSE VARIABLES ------------------------------------------------------

# __ Assigning experimental groups to individuals -------------------------
SuppA$Exp.Group <- NA
SuppA$Exp.Group[!is.na(SuppA$QC1x_1) | !is.na(SuppA$QC1x_2)] <- "Group A:\nGeneral\nproblem\nreminder"
SuppA$Exp.Group[!is.na(SuppA$QC2x_1) | !is.na(SuppA$QC2x_2)] <- "Group B:\nEconomic\nproblem"
SuppA$Exp.Group[!is.na(SuppA$QC3x_1) | !is.na(SuppA$QC3x_2)] <- "Group C:\nDeservingness\nproblem"
SuppA$Exp.Group[!is.na(SuppA$QC4x_1) | !is.na(SuppA$QC4x_2)] <- "Group D:\nEconomic &\ndeservingness\nproblem"
SuppA$Exp.Group[!is.na(SuppA$QC5x_1) | !is.na(SuppA$QC5x_2)] <- "Group E:\nControl\nGroup"
SuppA$Exp.Group <- factor(SuppA$Exp.Group)





# __ Merging responses of all experimental groups into a single column -----
# ____ SICKNESS benefits (variable name: QC*x_1) ---------------------------
SuppA$Sickness <- rowSums(SuppA[c(3, 5, 7, 9, 11)], na.rm = TRUE)
SuppA$Sickness[SuppA$Sickness == 0] <- NA

# Recoding "Don't know" into NA
SuppA$Sickness[SuppA$Sickness == 8] <- NA

# Adjusting range of the variable to <0;1>
SuppA$Sickness <- (SuppA$Sickness-1)/6





# ____ UNEMPLOYMENT benefits (variable name: QC*x_2) -----------------------
SuppA$Unemployment <- rowSums(SuppA[c(4, 6, 8, 10, 12)], na.rm = TRUE)
SuppA$Unemployment[SuppA$Unemployment == 0] <- NA

# Recoding "Don't know" into NA
SuppA$Unemployment[SuppA$Unemployment == 8] <- NA

# Adjusting range of the variable to <0;1>
SuppA$Unemployment <- (SuppA$Unemployment-1)/6


# Getting rid of the individual variables per experimental groups
SuppA[c(3:12)] <- NULL





# IDEOLOGY VARIABLES (Used for second part of analysis) ---------------------
# __ Left-right: As continuous  ---------------------------------------------
# Renaming variable
names(SuppA)[names(SuppA) == "Q34"] <- "Left.Right"     

# Recoding "Don't know" into NA
SuppA$Left.Right[SuppA$Left.Right == 12] <- NA          

# Setting the type
SuppA$Left.Right <- as.numeric(SuppA$Left.Right)

# Making the range from <0;1>
SuppA$Left.Right <- SuppA$Left.Right-1






# CONTROL VARIABLES (Used for robustness/randomization checks) --------------

# __ Gender  ----------------------------------------------------------------
# Renaming variable
names(SuppA)[names(SuppA) == "kjonn"] <- "Gender"

# Recoding values
SuppA$Gender[SuppA$Gender == 1] <- "Male"
SuppA$Gender[SuppA$Gender == 2] <- "Female"
SuppA$Gender <- factor(SuppA$Gender, levels = c("Male", "Female"))





# __ Age  -----------------------------------------
# Renaming variable
names(SuppA)[names(SuppA) == "alder"] <- "Age"

# Recoding "0" into NA
SuppA$Age[SuppA$Age == 0] <- NA  

# Adjusting format
SuppA$Age <- as.numeric(SuppA$Age)





# __ Education -----------------------------------------
# Renaming variable
names(SuppA)[names(SuppA) == "utdanning"] <- "Education"

# Recoding values
SuppA$Education[SuppA$Education == 1] <- "Primary school"
SuppA$Education[SuppA$Education == 2] <- "Higher education"
SuppA$Education[SuppA$Education == 3] <- "Vocational education"
SuppA$Education[SuppA$Education == 4] <- "University (up to 4 years)"
SuppA$Education[SuppA$Education == 5] <- "University (4+ years)"
SuppA$Education <- factor(SuppA$Education,
                          levels = c("Primary school",
                                     "Higher education",
                                     "Vocational education",
                                     "University (up to 4 years)",
                                     "University (4+ years)"))





# __ Income (Personal gross) -----------------------------------------
# CATEOGIRES: 1 - Below 200.000 kroner 
#             2 - 200.000 - 299.999 kroner
#             3 - 300.000 - 399.000 kroner
#             4 - 400.000 - 499.999 kroner
#             5 - 500.000 - 599.999 kroner
#             6 - 600.000 - 699.999 kroner
#             7 - 700-000 - 799.999 kroner
#             8 - 800.000 - 999.999 kroner
#             9 - 1.000.000 kroner or more
#             10- Decline to answer

# Renaming variable
names(SuppA)[names(SuppA) == "persinnt12"] <- "Income"

# Recoding "Decline to answer" into NA
SuppA$Income[SuppA$Income == 10] <- NA

# Adjusting format
SuppA$Income <- as.numeric(SuppA$Income)





# __ Weights -----------------------------------------
# Renaming variable
names(SuppA)[names(SuppA) == "vekt"] <- "Weight"

# Adjusting format
SuppA$Weight <- as.numeric(SuppA$Weight)










# FIGURE 1: DISTRIBUTION OF RESPONSES ---------------------------------------

# Creating a dataframe with experiment participants only
# (to make it easier to calculate group totals)
SuppA.summary <- subset(SuppA, !is.na(Exp.Group))

# Selecting only relevant variables
SuppA.summary <- SuppA.summary[c(1, 10:11)]

# Back-transforming the response variables
# (to get values as they were used in the survey)
SuppA.summary$Sickness <- (SuppA.summary$Sickness*6)+1
SuppA.summary$Unemployment <- (SuppA.summary$Unemployment*6)+1

# Getting distribution of individual responses
SuppA.summary <- 
  merge(
    data.frame(table(SuppA.summary$Sickness, useNA = "ifany")),
    data.frame(table(SuppA.summary$Unemployment, useNA = "ifany")),
    by = "Var1"
  )

# Changing names
names(SuppA.summary) <- c("survey.reply", "sickness", "unemployment")

# Adding shares
SuppA.summary$sickness.perc <- SuppA.summary$sickness/sum(SuppA.summary$sickness)
SuppA.summary$unemployment.perc <- SuppA.summary$unemployment/sum(SuppA.summary$unemployment)

# Creating labels
SuppA.summary$sickness.label.n <- paste("(", "italic('N')==", SuppA.summary$sickness, ")", sep = "")
SuppA.summary$unemployment.label.n <- paste("(", "italic('N')==", SuppA.summary$unemployment, ")", sep = "")

SuppA.summary$sickness.label.perc <- paste(round(SuppA.summary$sickness.perc*100, digits = 1), "%", sep="")
SuppA.summary$unemployment.label.perc <- paste0(round(SuppA.summary$unemployment.perc*100, digits = 1), "%", sep="")

# Minor adjustment (paste() command was erasing zero after decimal points)
SuppA.summary[1, "unemployment.label.perc"] <- "11.0%"



# Figure
figure <- 
  ggarrange(
    # Sickness benefits
    ggplot(data = na.omit(SuppA.summary), aes(x = survey.reply, y = sickness.perc)) +
      geom_hline(yintercept = c(0.1, 0.2), color = "grey80") +
      geom_bar(stat = "identity", width = 0.85, alpha = 0.5, colour = "black", linewidth = 0.5) +
      geom_text(aes(label = sickness.label.perc), vjust = -2, size = 3) +
      geom_text(aes(label = sickness.label.n), vjust = -0.15, size = 3, parse = TRUE) +
      scale_y_continuous(breaks = seq(0, 0.3, 0.02),
                         labels = c("0%", "", "", "", "", "10%", "", "", "", "", "20%", "", "", "", "", "30%")) +
      scale_x_discrete(limits = factor(1:7),
                       labels = c("1\nNorway will not\nbe able to afford\nthe present level", 
                                  "2\n\n\n", "3\n\n\n", "4\n\n\n", "5\n\n\n", "6\n\n\n", 
                                  "7\nNorway will be able\nto afford to increase\nthe present level")) +
      coord_cartesian(xlim = c(0.5, 7.5), ylim = c(0, 0.36), expand = FALSE, clip = "off") +
      labs(title = "Think ten years ahead in time. For each of the following social security systems\nand public services, where would you place yourself ...",
           subtitle = "   ",
           x = NULL, y = NULL) +
      geom_rect(xmin = 0.5, xmax = 7.5, ymin = 0.3, ymax = 0.36, colour = "black", fill = "grey80") +
      annotate("text", x = 4, y = 0.33, hjust = 0.5, vjust = 0.00, size = 5, label = "Sickness benefits", fontface = 2) +
      annotate("text", x = 4, y = 0.33, hjust = 0.5, vjust = 1.75, size = 3, label = paste0("[Std. Dev. = ", round(sd(na.omit(SuppA$Sickness)*6+1), 2),"]")) +
      theme(plot.title = element_text(size = 12, vjust = 0.5, hjust = 0.5),
            plot.subtitle = element_text(size = 5, vjust = 0.5, hjust = 0.5, color = "transparent"),
            strip.background = element_rect(colour = "black"),
            panel.grid.major.x = element_blank(), 
            panel.grid.major.y = element_line(colour = "gray90"),
            panel.grid.minor = element_blank(),
            panel.background = element_rect(fill = "transparent"),
            panel.border = element_rect(colour = "black", fill = "transparent"),
            axis.text.x = element_blank(),
            axis.ticks = element_blank(),
            plot.margin=unit(c(0.5, 0.5, 0.25, 0.5),"cm")),
    
    # Unemployment benefits
    ggplot(data = na.omit(SuppA.summary), aes(x = survey.reply, y = unemployment.perc)) +
      geom_hline(yintercept = c(0.1, 0.2), color = "grey80") +
      geom_bar(stat = "identity", width = 0.85, alpha = 0.5, colour = "black", linewidth = 0.5) +
      geom_text(aes(label = unemployment.label.perc), vjust = -2, size = 3) +
      geom_text(aes(label = unemployment.label.n), vjust = -0.15, size = 3, parse = TRUE) +
      scale_y_continuous(breaks = seq(0, 0.3, 0.02),
                         labels = c("0%", "", "", "", "", "10%", "", "", "", "", "20%", "", "", "", "", "30%")) +
      scale_x_discrete(limits = factor(1:7),
                       labels = c("1\nNorway will not\nbe able to afford\nthe present level", 
                                  "2\n\n\n", "3\n\n\n", "4\n\n\n", "5\n\n\n", "6\n\n\n", 
                                  "7\nNorway will be able\nto afford to increase\nthe present level")) +
      coord_cartesian(xlim = c(0.5, 7.5), ylim = c(0, 0.36), expand = FALSE, clip = "off") +
      labs(x = NULL, y = NULL) +
      geom_rect(xmin = 0.5, xmax = 7.5, ymin = 0.3, ymax = 0.36, colour = "black", fill = "grey80") +
      annotate("text", x = 4, y = 0.33, hjust = 0.5, vjust = 0.00, size = 5, label = "Unemployment benefits", fontface = 2) +
        annotate("text", x = 4, y = 0.33, hjust = 0.5, vjust = 1.75, size = 3, label = paste0("[Std. Dev. = ", round(sd(na.omit(SuppA$Unemployment)*6+1), 2),"]")) +
      theme(strip.background = element_rect(colour = "black"),
            panel.grid.major.x = element_blank(), 
            panel.grid.major.y = element_line(colour = "gray90"), 
            panel.grid.minor = element_blank(),
            panel.background = element_rect(fill = "transparent"),
            panel.border = element_rect(colour = "black", fill = "transparent"),
            axis.text.x = element_text(vjust = 0.5, hjust = 0.5),
            axis.ticks = element_blank(),
            plot.margin=unit(c(0.25, 0.5, 0.5, 0.5),"cm")),
    
    # Specs
    nrow = 2, heights = c(10, 10))


# Output
ggsave(figure, filename = "04-figures and models/Fig1 - histogram answers.png", width = 20, height = 17.5, units = "cm", dpi = 600)


# Cleanup
rm(figure, SuppA.summary)









# ANALYSIS: -----------------------------------------------------------------
# __ AVERAGE TREATMENT EFFECT -----------------------------------------------

# Setting 'Group E: Control Group' as a reference category
SuppA$Exp.Group <- factor(SuppA$Exp.Group, 
                          levels = c("Group E:\nControl\nGroup", 
                                     "Group A:\nGeneral\nproblem\nreminder", 
                                     "Group B:\nEconomic\nproblem", 
                                     "Group C:\nDeservingness\nproblem", 
                                     "Group D:\nEconomic &\ndeservingness\nproblem"))





# __ MODELS: Average treatment effects ------------------------------------
# Sickness benefits
summary(
  lm(Sickness ~ Exp.Group, 
     data = SuppA, 
     weights = Weight))


# Unemployment benefits
summary(
  lm(Unemployment ~ Exp.Group, 
     data = SuppA, 
     weights = Weight))





# ____ FIGURE 2: Forest plots ----------------------------------------------

ggpubr::ggarrange(
  
  # Strips
  grobTree(rectGrob(gp = gpar(fill = "grey90")), textGrob("Sickness benefits", hjust = 0.5,  gp = gpar(col = "black", cex = 1))),
  NULL,
  grobTree(rectGrob(gp = gpar(fill = "grey90")), textGrob("Unemployment benefits", hjust = 0.5,  gp = gpar(col = "black", cex = 1))),
  
  #Sickness benefits
  plot_model(lm(Sickness ~ Exp.Group, data = SuppA, weights = Weight),
             title = "",
             colors = "black",
             vline.color = "grey65") +
    scale_x_discrete(labels = rev(c("Group A:\nGeneral\nproblem\nreminder", 
                                    "Group B:\nEconomic\nproblem", 
                                    "Group C:\nDeservingness\nproblem", 
                                    "Group D:\nEconomic &\ndeservingness\nproblem")), 
                     expand = expansion(add = 0.3)) +
    scale_y_continuous(breaks = seq(from = -0.2, to = 0.2, by = 0.05), labels = c("-0.2", "", "-0.1", "", "0.0", "", "0.1", "", "0.2"),
                       limits = c(-0.2, 0.2), expand = expansion(add = 0)) +
    geom_hline(yintercept = c(-0.2, 0.2)) +
    geom_vline(xintercept = c(0.7, 4.3)) +
    labs(title = "Reference:\nControl group (E)") +
    theme(plot.title = element_text(size = 8, hjust = 0.5), 
          axis.text = element_text(size = 8, colour = "black"), axis.title = element_text(size = 9, colour = "black"),
          panel.grid.major.x = element_line(colour = "gray90"),
          panel.grid.major.y = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_rect(fill = "transparent", colour = "black"),
          plot.background = element_rect(fill = "transparent", colour = "black"),
          axis.ticks = element_blank(),
          plot.margin = unit(c(0.5, 0.5, 0.5, 0.5),"cm")),
  NULL,
  
  #Unemployment benefits
  plot_model(lm(Unemployment ~ Exp.Group, data = SuppA, weights = Weight),
             title = "",
             colors = "black",
             vline.color = "grey65") +
    scale_x_discrete(labels = rev(c("Group A:\nGeneral\nproblem\nreminder", 
                                    "Group B:\nEconomic\nproblem", 
                                    "Group C:\nDeservingness\nproblem", 
                                    "Group D:\nEconomic &\ndeservingness\nproblem")), 
                     expand = expansion(add = 0.3)) +
    scale_y_continuous(breaks = seq(from = -0.2, to = 0.2, by = 0.05), labels = c("-0.2", "", "-0.1", "", "0.0", "", "0.1", "", "0.2"),
                       limits = c(-0.2, 0.2), expand = expansion(add = 0)) +
    geom_hline(yintercept = c(-0.2, 0.2)) +
    geom_vline(xintercept = c(0.7, 4.3)) +
    labs(title = "Reference:\nControl group (E)") +
    theme(plot.title = element_text(size = 8, hjust = 0.5), 
          axis.text = element_text(size = 8, colour = "black"), axis.title = element_text(size = 9, colour = "black"),
          panel.grid.major.x = element_line(colour = "gray90"),
          panel.grid.major.y = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_rect(fill = "transparent", colour = "black"),
          plot.background = element_rect(fill = "transparent", colour = "black"),
          axis.ticks = element_blank(), 
          plot.margin = unit(c(0.5, 0.5, 0.5, 0.5),"cm")),
  
  ncol = 3, widths  = c(5, 0.25, 5),
  nrow = 2, heights = c(0.75, 10)) +
  theme(plot.margin = unit(c(0.1, 0.1, 0.1, 0.1),"cm"),
        panel.background = element_rect(fill = "white", colour = "white"),
        plot.background = element_rect(fill = "white", colour = "white"))


# Output
ggsave(filename = "04-figures and models/Fig2 - forest plots.png", width = 17.5, height = 12.5, units = "cm", dpi = 600)






# INTERACTION: AVERAGE EFFECT * IDEOLOGY --------------------------------------------------

# Setting 'Group E: Control Group' as a reference category
SuppA$Exp.Group <- factor(SuppA$Exp.Group, 
                          levels = c("Group E:\nControl\nGroup", 
                                     "Group A:\nGeneral\nproblem\nreminder", 
                                     "Group B:\nEconomic\nproblem", 
                                     "Group C:\nDeservingness\nproblem", 
                                     "Group D:\nEconomic &\ndeservingness\nproblem"))


# __ MODELS: Average treatment effects * Left-right ideology ------------------------------
# Sickness benefits
summary(
  lm(Sickness ~ Exp.Group * Left.Right, 
     data = SuppA, 
     weights = Weight))



# Unemployment benefits
summary(
  lm(Unemployment ~ Exp.Group * Left.Right, 
     data = SuppA, 
     weights = Weight))




# ____ TABLE 2: Regression output -----------------------------------
m1 <- lm(Sickness ~ Exp.Group, data = SuppA, weights = Weight)
m2 <- lm(Unemployment ~ Exp.Group, data = SuppA,  weights = Weight)
m3 <- lm(Sickness ~ Exp.Group * Left.Right, data = SuppA, weights = Weight)
m4 <- lm(Unemployment ~ Exp.Group * Left.Right, data = SuppA,  weights = Weight)

stargazer(#Main models: average treatment effect
          m1, m2,
          
          #Interactions with left-right ideology
          m3, m4,
          
          #Specifications
          type = "html", out = "04-figures and models/Tab2 - Regression models.html",
          intercept.top = TRUE, intercept.bottom = FALSE, no.space = TRUE,
          dep.var.labels = c("Sickness benefits",
                             "Unemployment benefits",
                             "Sickness benefits",
                             "Unemployment benefits"),
          covariate.labels = c("Constant",
                               "Group A: General problem reminder",
                               "Group B: Economic problem",
                               "Group C: Deservingness problem",
                               "Group D: Economic and deservingness problem",
                               "Left-Right",
                               "Group A * Left-Right",
                               "Group B * Left-Right",
                               "Group C * Left-Right",
                               "Group D * Left-Right"))


rm(m1, m2, m3, m4)




# ____ FIGURE 3: Marginal effects of treatments across ideological orientation ------------

# UPPER HALF OF THE FIGURE 3: Sickness benefits
# Generating the figure as a ggplot object
fig.sickness <- interplot(lm(Sickness ~ Exp.Group * Left.Right, data = SuppA, weights = Weight), var2 = "Left.Right", var1 = "Exp.Group") + theme(panel.spacing = unit(0.75, "lines"))

# Changing factor levels to improve strip labels
levels(fig.sickness$data$value)[levels(fig.sickness$data$value) == "Exp.GroupGroup A:\nGeneral\nproblem\nreminder"] <- "Group A:\nGeneral problem reminder"
levels(fig.sickness$data$value)[levels(fig.sickness$data$value) == "Exp.GroupGroup B:\nEconomic\nproblem"] <- "Group B:\nEconomic problem"
levels(fig.sickness$data$value)[levels(fig.sickness$data$value) == "Exp.GroupGroup C:\nDeservingness\nproblem"] <- "Group C:\nDeservingness problem"
levels(fig.sickness$data$value)[levels(fig.sickness$data$value) == "Exp.GroupGroup D:\nEconomic &\ndeservingness\nproblem"] <- "Group D:\nEconomic &\ndeservingness problem"

# Adding visual adjustments to the figure
fig.sickness <- fig.sickness +
  coord_cartesian(xlim = c(0, 10), ylim = c(-0.3, 0.3), expand = FALSE, clip = "off") +
  scale_y_continuous(breaks = seq(-0.3, 0.3, 0.1), labels = c("-0.3", "-0.2", "-0.1", "0.0", "0.1", "0.2", "0.3")) +
  scale_x_continuous(breaks = seq(0, 10, 1), labels = c("0", "\nLeft", "2", "", "4", "", "6", "", "8", "\nRight", "10")) +
  geom_hline(yintercept = 0, color = brewer.pal(n = 3, name = "Set1")[1], linetype = "longdash", alpha = 0.5) +
  geom_text(data = data.frame(value = c("Group A:\nGeneral problem reminder"), label = c("Reference: Control group (E)")), 
            aes(x = 0, y = 0, label = label), hjust = -0.025, vjust = 1.5, 
            size = 2.75, color = brewer.pal(n = 3, name = "Set1")[1], alpha = 0.5, 
            inherit.aes = FALSE) +
  theme(strip.background = element_rect(fill = "gray90", colour = "black"),
        panel.grid.major.x = element_line(colour = "gray90"), panel.grid.minor.x = element_blank(),
        panel.grid.major.y = element_line(colour = "gray90"), panel.grid.minor.y = element_line(colour = "gray90"),
        panel.background = element_rect(fill = "transparent", colour = "black"),
        panel.border = element_rect(fill = "transparent", colour = "black", size = 0.5),
        axis.text.x = element_blank(),
        axis.ticks = element_blank(),
        plot.margin = unit(c(0, -0.1, 0.15, 0.25),"cm"))

# Adding the right stripe to the figure
fig.sickness <- 
  egg::ggarrange(fig.sickness,
                 
                 # right hand stripe specs
                 ggplot(data.frame(a = c(0, 10), b = c(-0.2, 0.2)), aes(x = a, y = b)) +
                   coord_cartesian(xlim = c(0, 10), ylim = c(-0.2, 0.2), expand = FALSE, clip = "off") +
                   geom_rect(xmin = 0, xmax = 10, ymin = -0.2, ymax = 0.2, colour = "black", fill = "grey90", size = 0.5) +
                   annotate("text", x = 5, y = 0, hjust = 0.5, vjust = 0.35, angle = 270, size = 3.5, label = "Sickness benefits") +
                   theme(text = element_blank(), line = element_blank(), plot.margin = unit(c(0, 0.5, 0.15, 0),"cm")),
                 ncol = 2, widths = c(10, 0.5)
  ) 



# BOTTOM HALF OF THE FIGURE 3: Sickness benefits
# Generating the figure as a ggplot object
fig.unemployment <- interplot(lm(Unemployment ~ Exp.Group * Left.Right, data = SuppA, weights = Weight), var2 = "Left.Right", var1 = "Exp.Group") + theme(panel.spacing = unit(0.75, "lines"))

# Changing factor levels to improve strip labels
levels(fig.unemployment$data$value)[levels(fig.unemployment$data$value) == "Exp.GroupGroup A:\nGeneral\nproblem\nreminder"] <- "Group A:\nGeneral problem reminder"
levels(fig.unemployment$data$value)[levels(fig.unemployment$data$value) == "Exp.GroupGroup B:\nEconomic\nproblem"] <- "Group B:\nEconomic problem"
levels(fig.unemployment$data$value)[levels(fig.unemployment$data$value) == "Exp.GroupGroup C:\nDeservingness\nproblem"] <- "Group C:\nDeservingness problem"
levels(fig.unemployment$data$value)[levels(fig.unemployment$data$value) == "Exp.GroupGroup D:\nEconomic &\ndeservingness\nproblem"] <- "Group D:\nEconomic &\ndeservingness problem"

# Adding visual adjustments to the figure
fig.unemployment <- fig.unemployment +
  coord_cartesian(xlim = c(0, 10), ylim = c(-0.3, 0.3), expand = FALSE, clip = "off") +
  scale_y_continuous(breaks = seq(-0.3, 0.3, 0.1), labels = c("-0.3", "-0.2", "-0.1", "0.0", "0.1", "0.2", "0.3")) +
  scale_x_continuous(breaks = seq(0, 10, 1), labels = c("0", "\nLeft", "2", "", "4", "", "6", "", "8", "\nRight", "10")) +
  geom_hline(yintercept = 0, color = brewer.pal(n = 3, name = "Set1")[1], linetype = "longdash", alpha = 0.5) +
  geom_text(data = data.frame(value = c("Group A:\nGeneral problem reminder"), label = c("Reference: Control group (E)")), 
            aes(x = 0, y = 0, label = label), hjust = -0.025, vjust = 1.5, 
            size = 2.75, color = brewer.pal(n = 3, name = "Set1")[1], alpha = 0.5, 
            inherit.aes = FALSE) +
  theme(strip.background = element_blank(), strip.text = element_blank(),
        panel.grid.major.x = element_line(colour = "gray90"), panel.grid.minor.x = element_blank(),
        panel.grid.major.y = element_line(colour = "gray90"), panel.grid.minor.y = element_line(colour = "gray90"),
        panel.background = element_rect(fill = "transparent", colour = "black"),
        panel.border = element_rect(fill = "transparent", colour = "black"),
        axis.ticks = element_blank(),
        plot.margin = unit(c(0, -0.1, 0.15, 0.25),"cm"))


# Adding the right stripe to the figure
fig.unemployment <- 
  egg::ggarrange(fig.unemployment,
                 
                 # right hand stripe specs
                 ggplot(data.frame(a = c(0, 10), b = c(-0.2, 0.2)), aes(x = a, y = b)) +
                   coord_cartesian(xlim = c(0, 10), ylim = c(-0.2, 0.2), expand = FALSE, clip = "off") +
                   geom_rect(xmin = 0, xmax = 10, ymin = -0.2, ymax = 0.2, colour = "black", fill = "grey90", size = 0.5) +
                   annotate("text", x = 5, y = 0, hjust = 0.5, vjust = 0.35, angle = 270, size = 3.5, label = "Unemployment benefits") +
                   theme(text = element_blank(), line = element_blank(), plot.margin = unit(c(0, 0.5, 0, 0),"cm")),
                 ncol = 2, widths = c(10, 0.5)
  )




# PUTTING THE TWO HALVES TOGETHER
annotate_figure(ggpubr::ggarrange(fig.sickness, fig.unemployment, nrow = 2, heights = c(11, 10)),
                bottom = "Political ideology",
                left = "Marginal effect") + theme(plot.margin = unit(c(0.25, 0, 0.25, 0.25),"cm")) +
  theme(plot.margin = unit(c(0.1, 0.1, 0.1, 0.1),"cm"),
        panel.background = element_rect(fill = "white", colour = "white"),
        plot.background = element_rect(fill = "white", colour = "white"))


# Saving the output
ggsave(filename = "04-figures and models/Fig3 - Marginal plot - effect vs ideology.png", width = 25, height = 12.5, units = "cm", dpi = 600)

rm(fig.sickness, fig.unemployment)









# APPENDIX MATERIAL -----------------------------------------------------------------------

# __ TABLE A2: NUMBER OF PARTICIPANTS PER EXPERIMENTAL GROUPS -----------------------------

# Creating a data.frame with experiment participants only
# (to make it easier to calculate group totals)
SuppA.summary <- subset(SuppA, !is.na(Exp.Group))


# Creating identification dummy
SuppA.summary$Exp.sickness <- ifelse(is.na(SuppA.summary$Sickness), 0, 1)
SuppA.summary$Exp.Unemployment <- ifelse(is.na(SuppA.summary$Unemployment), 0, 1)


# Table with group total
summary.groups <- 
  cbind(
    subset(data.frame(table(SuppA.summary$Exp.Group, SuppA.summary$Exp.sickness)), Var2 == 1),
    subset(data.frame(table(SuppA.summary$Exp.Group, SuppA.summary$Exp.Unemployment)), Var2 == 1)
  )[c(1, 3, 6)]

names(summary.groups) <- c("Group", "Sickness benefits", "Unemployment benefits")


# Transposing the table
SuppA.summary <- dcast(reshape2::melt(summary.groups, id.vars = "Group"), variable ~ Group)


# Adjusting the column names
names(SuppA.summary) <- c("Experiment - part", "Group E: Control Group", "Group A: General problem reminder", "Group B: Economic problem", "Group C: Deservingness problem", "Group D: Economic and deservingness problem")



# Table output
write.csv(SuppA.summary, file = "04-figures and models/TabA2 - Group totals.csv")

# Cleanup
rm(SuppA.summary, summary.groups)






# __ MODELS with controls: Average treatment effects  ---------------------

# Setting 'Group E: Control Group' as a reference category
SuppA$Exp.Group <- factor(SuppA$Exp.Group, 
                          levels = c("Group E:\nControl\nGroup",
                                     "Group B:\nEconomic\nproblem", 
                                     "Group C:\nDeservingness\nproblem", 
                                     "Group D:\nEconomic &\ndeservingness\nproblem",
                                     "Group A:\nGeneral\nproblem\nreminder"))



# Models: Average treatment effects with controls
# Sickness benefits
m1 <- lm(Sickness ~ Exp.Group 
                          + as.factor(Gender) + Age + as.factor(Education) + as.factor(Income), 
                          data = SuppA, weights = Weight)

summary(m1)


# Unemployment benefits
m2 <- lm(Unemployment ~ Exp.Group
                              + as.factor(Gender) + Age + as.factor(Education) + as.factor(Income), 
                              data = SuppA, weights = Weight)

summary(m2)



# Models: Treatment effects along left-right dimension with controls
# Sickness benefits
m3 <- lm(Sickness ~ Exp.Group * Left.Right
                          + as.factor(Gender) + Age + as.factor(Education) + as.factor(Income), 
                          data = SuppA, weights = Weight)

summary(m3)


# Unemployment benefits
m4 <- lm(Unemployment ~ Exp.Group  * Left.Right
                              + as.factor(Gender) + Age + as.factor(Education) + as.factor(Income), 
                              data = SuppA, weights = Weight)

summary(m4)


# ____ TABLE A3: Regression output -----------------------------------
stargazer(m1, m2,
          m3, m4,
          type = "html", out = "04-figures and models/TabA3 - Models with controls.html",
          intercept.top = TRUE, intercept.bottom = FALSE, no.space = TRUE,
          dep.var.labels = c("Sickness benefits",
                             "Unemployment benefits",
                             "Sickness benefits",
                             "Unemployment benefits"),
          covariate.labels = c("Constant",
                               "Group B: Economic problem",
                               "Group C: Deservingness problem",
                               "Group D: Economic and deservingness problem",
                               "Group A: Problem reminder",
                               "Left-Right",
                               "Female",
                               "Age",
                               "Higher education",
                               "Vocational education",
                               "University (up to 4 years)",
                               "University (4+ years)",
                               "Income: 200 000-299 999 kr",
                               "300 000-399 000 kr",
                               "400 000-499 999 kr",
                               "500 000-599 999 kr",
                               "600 000-699 999 kr",
                               "700 000-799 999 kr",
                               "800 000-999 999 kr",
                               "1 000 000 kr or more",
                               "Group B * Left-Right", 
                               "Group C * Left-Right", 
                               "Group D * Left-Right",
                               "Group A * Left-Right"))


rm(m1, m2,
   m3, m4)






# __ TABLE A4: COMPARISON OF EXPERIMENTAL GROUPS ---------------------------------------
# I.e., checking whether randomization worked as expected


# ____ Sickness benefits ------------------------------------------------
# GENDER
# Subsetting the participants
SuppA.summary <- subset(SuppA, !is.na(Exp.Group) & !is.na(Sickness))

# Exporting frequencies
Gender <- data.frame(table(SuppA.summary$Exp.Group, SuppA.summary$Gender, useNA = "always"))
Gender <- subset(Gender, !is.na(Gender$Var1))

# Adding totals
Gender <- merge(Gender, aggregate(Gender$Freq, by = list(Category = Gender$Var1), FUN = sum), by.x = "Var1", by.y = "Category")

# Calculating shares
Gender$share <- round((Gender$Freq/Gender$x)*100, digits = 1)

# Changing names
names(Gender) <- c("Group", "Category", "N", "Group total", "%")

# Adding variable name and changing the order
Gender$Variable <- c("Gender")
Gender <- Gender[c(1, 6, 2:5)]



# AGE
# Subsetting the participants
SuppA.summary <- subset(SuppA, !is.na(Exp.Group) & !is.na(Sickness))

# Exporting averages
Age <- aggregate(SuppA.summary$Age, by = list(Category = SuppA.summary$Exp.Group), FUN = mean)

# Putting data.frame into a consistent shape
Age <- data.frame("Group"       = Age$Category,
                  "Variable"    = c("Age"),
                  "Category"    = c(NA),
                  "N"           = c(NA),
                  "Group total" = c(NA),
                  "%"           = Age$x)

# Changing names
names(Age) <- c("Group", "Variable", "Category", "N", "Group total", "%")



# EDUCATION
# Subsetting the participants
SuppA.summary <- subset(SuppA, !is.na(Exp.Group) & !is.na(Sickness))

# Exporting frequencies
Education <- data.frame(table(SuppA.summary$Exp.Group, SuppA.summary$Education, useNA = "always"))
Education <- subset(Education, !is.na(Education$Var1))

# Adding totals
Education <- merge(Education, aggregate(Education$Freq, by = list(Category = Education$Var1), FUN = sum), by.x = "Var1", by.y = "Category")

# Calculating shares
Education$share <- round((Education$Freq/Education$x)*100, digits = 1)

# Changing names
names(Education) <- c("Group", "Category", "N", "Group total", "%")

# Adding variable name and changing the order
Education$Variable <- c("Education")
Education <- Education[c(1, 6, 2:5)]



# INCOME
# Subsetting the participants
SuppA.summary <- subset(SuppA, !is.na(Exp.Group) & !is.na(Sickness))

# Exporting frequencies
Income <- data.frame(table(SuppA.summary$Exp.Group, as.factor(SuppA.summary$Income), useNA = "always"))
Income <- subset(Income, !is.na(Income$Var1))

# Adding totals
Income <- merge(Income, aggregate(Income$Freq, by = list(Category = Income$Var1), FUN = sum), by.x = "Var1", by.y = "Category")

# Calculating shares
Income$share <- round((Income$Freq/Income$x)*100, digits = 1)

# Changing names
names(Income) <- c("Group", "Category", "N", "Group total", "%")

# Adding variable name and changing the order
Income$Variable <- c("Income")
Income <- Income[c(1, 6, 2:5)]



# Merging individual summaries together and creating a wide table
Summary.table <- rbind(Gender, Age, Education, Income)
rm(Gender, Age, Education, Income, SuppA.summary)

# Creating a wide dataset
Summary.table <-
  cbind(
    subset(Summary.table, Group == "Group A:\nGeneral\nproblem\nreminder")[order(subset(Summary.table, Group == "Group A:\nGeneral\nproblem\nreminder")$Variable, subset(Summary.table, Group == "Group A:\nGeneral\nproblem\nreminder")$Category), ],
    subset(Summary.table, Group == "Group B:\nEconomic\nproblem")[order(subset(Summary.table, Group == "Group B:\nEconomic\nproblem")$Variable, subset(Summary.table, Group == "Group B:\nEconomic\nproblem")$Category), ],
    subset(Summary.table, Group == "Group C:\nDeservingness\nproblem")[order(subset(Summary.table, Group == "Group C:\nDeservingness\nproblem")$Variable, subset(Summary.table, Group == "Group C:\nDeservingness\nproblem")$Category), ],
    subset(Summary.table, Group == "Group D:\nEconomic &\ndeservingness\nproblem")[order(subset(Summary.table, Group == "Group D:\nEconomic &\ndeservingness\nproblem")$Variable, subset(Summary.table, Group == "Group D:\nEconomic &\ndeservingness\nproblem")$Category), ],
    subset(Summary.table, Group == "Group E:\nControl\nGroup")[order(subset(Summary.table, Group == "Group E:\nControl\nGroup")$Variable, subset(Summary.table, Group == "Group E:\nControl\nGroup")$Category), ]
  )

# Leaving relevant information
Summary.table <- Summary.table[c(2:3, 4, 6, 10, 12, 16, 18, 22, 24, 28, 30)]
names(Summary.table) <- c("Variable", "Category", "Group A: N", "Group A: %", "Group B: N", "Group B: %", "Group C: N", "Group C: %", "Group D: N", "Group D: %", "Group E: N", "Group E: %")

# Rounding shares
Summary.table[c(4, 6, 8, 10, 12)] <- round(Summary.table[c(4, 6, 8, 10, 12)], digits = 1)

# Adding totals
Summary.table[21, ] <- c("Group total N", "", 
                         sum(subset(Summary.table, Variable == "Gender")[c( 3)]), "", 
                         sum(subset(Summary.table, Variable == "Gender")[c( 5)]), "",
                         sum(subset(Summary.table, Variable == "Gender")[c( 7)]), "",
                         sum(subset(Summary.table, Variable == "Gender")[c( 9)]), "",
                         sum(subset(Summary.table, Variable == "Gender")[c(11)]), "")

write.csv(Summary.table, file = "04-figures and models/TabA5a - Group descriptives - Sickness benefits.csv")
rm(Summary.table)




# ____ Unemployment benefits ------------------------------------------------
# GENDER
# Subsetting the participants
SuppA.summary <- subset(SuppA, !is.na(Exp.Group) & !is.na(Unemployment))

# Exporting frequencies
Gender <- data.frame(table(SuppA.summary$Exp.Group, SuppA.summary$Gender, useNA = "always"))
Gender <- subset(Gender, !is.na(Gender$Var1))

# Adding totals
Gender <- merge(Gender, aggregate(Gender$Freq, by = list(Category = Gender$Var1), FUN = sum), by.x = "Var1", by.y = "Category")

# Calculating shares
Gender$share <- round((Gender$Freq/Gender$x)*100, digits = 1)

# Changing names
names(Gender) <- c("Group", "Category", "N", "Group total", "%")

# Adding variable name and changing the order
Gender$Variable <- c("Gender")
Gender <- Gender[c(1, 6, 2:5)]



# AGE
# Subsetting the participants
SuppA.summary <- subset(SuppA, !is.na(Exp.Group) & !is.na(Unemployment))

# Exporting averages
Age <- aggregate(SuppA.summary$Age, by = list(Category = SuppA.summary$Exp.Group), FUN = mean)

# Putting data.frame into a consistent shape
Age <- data.frame("Group"       = Age$Category,
                  "Variable"    = c("Age"),
                  "Category"    = c(NA),
                  "N"           = c(NA),
                  "Group total" = c(NA),
                  "%"           = Age$x)

# Changing names
names(Age) <- c("Group", "Variable", "Category", "N", "Group total", "%")



# EDUCATION
# Subsetting the participants
SuppA.summary <- subset(SuppA, !is.na(Exp.Group) & !is.na(Unemployment))

# Exporting frequencies
Education <- data.frame(table(SuppA.summary$Exp.Group, SuppA.summary$Education, useNA = "always"))
Education <- subset(Education, !is.na(Education$Var1))

# Adding totals
Education <- merge(Education, aggregate(Education$Freq, by = list(Category = Education$Var1), FUN = sum), by.x = "Var1", by.y = "Category")

# Calculating shares
Education$share <- round((Education$Freq/Education$x)*100, digits = 1)

# Changing names
names(Education) <- c("Group", "Category", "N", "Group total", "%")

# Adding variable name and changing the order
Education$Variable <- c("Education")
Education <- Education[c(1, 6, 2:5)]



# INCOME
# Subsetting the participants
SuppA.summary <- subset(SuppA, !is.na(Exp.Group) & !is.na(Unemployment))

# Exporting frequencies
Income <- data.frame(table(SuppA.summary$Exp.Group, as.factor(SuppA.summary$Income), useNA = "always"))
Income <- subset(Income, !is.na(Income$Var1))

# Adding totals
Income <- merge(Income, aggregate(Income$Freq, by = list(Category = Income$Var1), FUN = sum), by.x = "Var1", by.y = "Category")

# Calculating shares
Income$share <- round((Income$Freq/Income$x)*100, digits = 1)

# Changing names
names(Income) <- c("Group", "Category", "N", "Group total", "%")

# Adding variable name and changing the order
Income$Variable <- c("Income")
Income <- Income[c(1, 6, 2:5)]



# Merging individual summaries together and creating a wide table
Summary.table <- rbind(Gender, Age, Education, Income)
rm(Gender, Age, Education, Income, SuppA.summary)

# Creating a wide dataset
Summary.table <-
  cbind(
    subset(Summary.table, Group == "Group A:\nGeneral\nproblem\nreminder")[order(subset(Summary.table, Group == "Group A:\nGeneral\nproblem\nreminder")$Variable, subset(Summary.table, Group == "Group A:\nGeneral\nproblem\nreminder")$Category), ],
    subset(Summary.table, Group == "Group B:\nEconomic\nproblem")[order(subset(Summary.table, Group == "Group B:\nEconomic\nproblem")$Variable, subset(Summary.table, Group == "Group B:\nEconomic\nproblem")$Category), ],
    subset(Summary.table, Group == "Group C:\nDeservingness\nproblem")[order(subset(Summary.table, Group == "Group C:\nDeservingness\nproblem")$Variable, subset(Summary.table, Group == "Group C:\nDeservingness\nproblem")$Category), ],
    subset(Summary.table, Group == "Group D:\nEconomic &\ndeservingness\nproblem")[order(subset(Summary.table, Group == "Group D:\nEconomic &\ndeservingness\nproblem")$Variable, subset(Summary.table, Group == "Group D:\nEconomic &\ndeservingness\nproblem")$Category), ],
    subset(Summary.table, Group == "Group E:\nControl\nGroup")[order(subset(Summary.table, Group == "Group E:\nControl\nGroup")$Variable, subset(Summary.table, Group == "Group E:\nControl\nGroup")$Category), ]
  )

# Leaving relevant information
Summary.table <- Summary.table[c(2:3, 4, 6, 10, 12, 16, 18, 22, 24, 28, 30)]
names(Summary.table) <- c("Variable", "Category", "Group A: N", "Group A: %", "Group B: N", "Group B: %", "Group C: N", "Group C: %", "Group D: N", "Group D: %", "Group E: N", "Group E: %")

# Rounding shares
Summary.table[c(4, 6, 8, 10, 12)] <- round(Summary.table[c(4, 6, 8, 10, 12)], digits = 1)

# Adding totals
Summary.table[21, ] <- c("Group total N", "", 
                         sum(subset(Summary.table, Variable == "Gender")[c( 3)]), "", 
                         sum(subset(Summary.table, Variable == "Gender")[c( 5)]), "",
                         sum(subset(Summary.table, Variable == "Gender")[c( 7)]), "",
                         sum(subset(Summary.table, Variable == "Gender")[c( 9)]), "",
                         sum(subset(Summary.table, Variable == "Gender")[c(11)]), "")

write.csv(Summary.table, file = "04-figures and models/TabA5b - Group descriptives - Unemployment benefits.csv")
rm(Summary.table)






# __ ROBUSTNESS-CHECK: USING CATEGORICAL LEFT-RIGHT SCALE ---------------------------------------

# Creating a categorical left-right variable
SuppA$Left.Right.Cat <- NA
SuppA$Left.Right.Cat[SuppA$Left.Right <= 2] <- "Left [0-2]"
SuppA$Left.Right.Cat[SuppA$Left.Right == 3 | SuppA$Left.Right == 4] <- "Centre Left [3-4]"
SuppA$Left.Right.Cat[SuppA$Left.Right == 5] <- "Center [5]"
SuppA$Left.Right.Cat[SuppA$Left.Right == 6 | SuppA$Left.Right == 7] <- "Centre Right [6-7]"
SuppA$Left.Right.Cat[SuppA$Left.Right >= 8] <- "Right [8-10]"


SuppA$Left.Right.Cat <- factor(SuppA$Left.Right.Cat, levels = c("Left [0-2]",
                                                                "Centre Left [3-4]",
                                                                "Center [5]",
                                                                "Centre Right [6-7]",
                                                                "Right [8-10]"))



# Setting 'Group E: Control Group' as a reference category
SuppA$Exp.Group <- factor(SuppA$Exp.Group, 
                          levels = c("Group E:\nControl\nGroup", 
                                     "Group A:\nGeneral\nproblem\nreminder", 
                                     "Group B:\nEconomic\nproblem", 
                                     "Group C:\nDeservingness\nproblem", 
                                     "Group D:\nEconomic &\ndeservingness\nproblem"))


# __ MODELS: Average treatment effects * Left-right ideology ------------------------------
# Sickness benefits
summary(
  lm(Sickness ~ Exp.Group * Left.Right.Cat, 
     data = SuppA, 
     weights = Weight))



# Unemployment benefits
summary(
  lm(Unemployment ~ Exp.Group * Left.Right.Cat, 
     data = SuppA, 
     weights = Weight))




# ____ TABLE A4: Regression output -----------------------------------
m1 <- lm(Sickness ~ Exp.Group * Left.Right.Cat, data = SuppA, weights = Weight)
m2 <- lm(Unemployment ~ Exp.Group * Left.Right.Cat, data = SuppA,  weights = Weight)

stargazer(m1, m2,
          # Specs
          type = "html", out = "04-figures and models/TabA4 - Alternative left-right.html",
          intercept.top = TRUE, intercept.bottom = FALSE, no.space = TRUE, single.row = TRUE,
          dep.var.labels = c("Sickness benefits",
                             "Unemployment benefits"),
          covariate.labels = c("Constant",
                               "Group A: General problem reminder",
                               "Group B: Economic problem",
                               "Group C: Deservingness problem",
                               "Group D: Economic and deservingness problem",
                               "Left-Right: Centre left [3-4]",
                               "Left-Right: Centre [5]",
                               "Left-Right: Centre right [6-7]",
                               "Left-Right: Right [8-10]",
                               "Group A * Centre left [3-4]",
                               "Group B * Centre left [3-4]",
                               "Group C * Centre left [3-4]",
                               "Group D * Centre left [3-4]",
                               "Group A * Centre [5]",
                               "Group B * Centre [5]",
                               "Group C * Centre [5]",
                               "Group D * Centre [5]",
                               "Group A * Centre right [6-7]",
                               "Group B * Centre right [6-7]",
                               "Group C * Centre right [6-7]",
                               "Group D * Centre right [6-7]",
                               "Group A * Right [8-10]",
                               "Group B * Right [8-10]",
                               "Group C * Right [8-10]",
                               "Group D * Right [8-10]"))

rm(m1, m2)





# ____ FIGURE A1: Marginal effects of treatments across ideological orientation ------------

ggpubr::ggarrange(
  
  # SICKNESS BENEFITS
  ggarrange(
    # Sickness benefits * Left ideology
    plot_model(lm(Sickness ~ Exp.Group, data = subset(SuppA, Left.Right.Cat == "Left [0-2]"), weights = Weight),
               type = "est",
               title = "",
               colors = "black",
               vline.color = brewer.pal(n = 9, name = "Paired")[5]) +
      geom_hline(yintercept = c(-0.30, -0.25, -0.20, -0.15, -0.10, -0.05, 0.05, 0.10, 0.15, 0.20, 0.25, 0.30), colour = "gray90") +
      geom_point(size = 2.5) +
      geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0, size = 0.7) +
      scale_x_discrete(limits = rev(c("Exp.GroupGroup A:\nGeneral\nproblem\nreminder", 
                                  "Exp.GroupGroup B:\nEconomic\nproblem", 
                                  "Exp.GroupGroup C:\nDeservingness\nproblem", 
                                  "Exp.GroupGroup D:\nEconomic &\ndeservingness\nproblem")),
                       labels = rev(c("Group A: General\nproblem reminder", 
                                      "Group B:\nEconomic problem", 
                                      "Group C:\nDeservingness problem", 
                                      "Group D: Economic &\ndeservingness problem"))) +
      scale_y_continuous(breaks = seq(from = -0.3, to = 0.3, by = 0.05), 
                         labels = c("-0.3", "", "-0.2", "", "-0.1", "", "0.0", "", "0.1", "", "0.2", "", "0.3")) +
      labs(x = NULL, y = NULL, title = NULL) +
      geom_rect(ymin = -0.3, ymax = 0.3, xmin = 4.5, xmax = 5.1, colour = "black", fill = "grey90") +
      annotate("text", x = 4.8, y = 0, hjust = 0.5, vjust = 0.5, label = "Sickness benefits", size = 11/.pt) +
      annotate("text", x = 3.5, y = -0.0125, hjust = 1, vjust = 0.5, label = "Reference:\nGroup E", size = 7/.pt, lineheight = 0.8, color = brewer.pal(n = 9, name = "Paired")[5]) +
      geom_hline(yintercept = c(-0.3, 0.3)) +
      geom_vline(xintercept = c(0.5, 4.5)) +
      coord_flip(ylim = c(-0.3, 0.3), xlim = c(0.5, 5.0), expand = FALSE, clip = "off") +
      theme(panel.grid = element_blank(),
            panel.background = element_rect(fill = "transparent", colour = "transparent"),
            plot.background = element_rect(fill = "transparent", colour = "transparent"),
            axis.ticks = element_blank(),
            axis.text.x = element_blank(),
            plot.margin = unit(c(0.5, 0.2, 0.2, 0.5),"cm")),
    
    # Sickness benefits * Left ideology
    plot_model(lm(Sickness ~ Exp.Group, data = subset(SuppA, Left.Right.Cat == "Centre Left [3-4]"), weights = Weight),
               type = "est",
               title = "",
               colors = "black",
               vline.color = brewer.pal(n = 9, name = "Paired")[5]) +
      geom_hline(yintercept = c(-0.30, -0.25, -0.20, -0.15, -0.10, -0.05, 0.05, 0.10, 0.15, 0.20, 0.25, 0.30), colour = "gray90") +
      geom_point(size = 2.5) +
      geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0, size = 0.7) +
      scale_x_discrete(limits = rev(c("Exp.GroupGroup A:\nGeneral\nproblem\nreminder", 
                                      "Exp.GroupGroup B:\nEconomic\nproblem", 
                                      "Exp.GroupGroup C:\nDeservingness\nproblem", 
                                      "Exp.GroupGroup D:\nEconomic &\ndeservingness\nproblem")),
                       labels = rev(c("Group A: General\nproblem reminder", 
                                      "Group B:\nEconomic problem", 
                                      "Group C:\nDeservingness problem", 
                                      "Group D: Economic &\ndeservingness problem"))) +
      scale_y_continuous(breaks = seq(from = -0.3, to = 0.3, by = 0.05), 
                         labels = c("-0.3", "", "-0.2", "", "-0.1", "", "0.0", "", "0.1", "", "0.2", "", "0.3")) +
      labs(x = NULL, y = NULL, title = NULL) +
      geom_hline(yintercept = c(-0.3, 0.3)) +
      geom_vline(xintercept = c(0.5, 4.5)) +
      coord_flip(ylim = c(-0.3, 0.3), xlim = c(0.5, 4.5), expand = FALSE, clip = "off") +
      theme(panel.grid = element_blank(),
            panel.background = element_rect(fill = "transparent", colour = "transparent"),
            plot.background = element_rect(fill = "transparent", colour = "transparent"),
            axis.ticks = element_blank(),
            axis.text.x = element_blank(),
            plot.margin = unit(c(0.5, 0.2, 0.2, 0.5),"cm")),
    
    
    # Sickness benefits * Center ideology
    plot_model(lm(Sickness ~ Exp.Group, data = subset(SuppA, Left.Right.Cat == "Center [5]"), weights = Weight),
               type = "est",
               title = "",
               colors = "black",
               vline.color = brewer.pal(n = 9, name = "Paired")[5]) +
      geom_hline(yintercept = c(-0.30, -0.25, -0.20, -0.15, -0.10, -0.05, 0.05, 0.10, 0.15, 0.20, 0.25, 0.30), colour = "gray90") +
      geom_point(size = 2.5) +
      geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0, size = 0.7) +
      scale_x_discrete(limits = rev(c("Exp.GroupGroup A:\nGeneral\nproblem\nreminder", 
                                      "Exp.GroupGroup B:\nEconomic\nproblem", 
                                      "Exp.GroupGroup C:\nDeservingness\nproblem", 
                                      "Exp.GroupGroup D:\nEconomic &\ndeservingness\nproblem")),
                       labels = rev(c("Group A: General\nproblem reminder", 
                                      "Group B:\nEconomic problem", 
                                      "Group C:\nDeservingness problem", 
                                      "Group D: Economic &\ndeservingness problem"))) +
      scale_y_continuous(breaks = seq(from = -0.3, to = 0.3, by = 0.05), 
                         labels = c("-0.3", "", "-0.2", "", "-0.1", "", "0.0", "", "0.1", "", "0.2", "", "0.3")) +
      labs(x = NULL, y = NULL, title = NULL) +
      geom_hline(yintercept = c(-0.3, 0.3)) +
      geom_vline(xintercept = c(0.5, 4.5)) +
      coord_flip(ylim = c(-0.3, 0.3), xlim = c(0.5, 4.5), expand = FALSE, clip = "off") +
      theme(panel.grid = element_blank(),
            panel.background = element_rect(fill = "transparent", colour = "transparent"),
            plot.background = element_rect(fill = "transparent", colour = "transparent"),
            axis.ticks = element_blank(),
            axis.text.x = element_blank(),
            plot.margin = unit(c(0.2, 0.2, 0.2, 0.5),"cm")),
    
    
    # Sickness benefits * Center ideology
    plot_model(lm(Sickness ~ Exp.Group, data = subset(SuppA, Left.Right.Cat == "Centre Right [6-7]"), weights = Weight),
               type = "est",
               title = "",
               colors = "black",
               vline.color = brewer.pal(n = 9, name = "Paired")[5]) +
      geom_hline(yintercept = c(-0.30, -0.25, -0.20, -0.15, -0.10, -0.05, 0.05, 0.10, 0.15, 0.20, 0.25, 0.30), colour = "gray90") +
      geom_point(size = 2.5) +
      geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0, size = 0.7) +
      scale_x_discrete(limits = rev(c("Exp.GroupGroup A:\nGeneral\nproblem\nreminder", 
                                      "Exp.GroupGroup B:\nEconomic\nproblem", 
                                      "Exp.GroupGroup C:\nDeservingness\nproblem", 
                                      "Exp.GroupGroup D:\nEconomic &\ndeservingness\nproblem")),
                       labels = rev(c("Group A: General\nproblem reminder", 
                                      "Group B:\nEconomic problem", 
                                      "Group C:\nDeservingness problem", 
                                      "Group D: Economic &\ndeservingness problem"))) +
      scale_y_continuous(breaks = seq(from = -0.3, to = 0.3, by = 0.05), 
                         labels = c("-0.3", "", "-0.2", "", "-0.1", "", "0.0", "", "0.1", "", "0.2", "", "0.3")) +
      labs(x = NULL, y = NULL, title = NULL) +
      geom_hline(yintercept = c(-0.3, 0.3)) +
      geom_vline(xintercept = c(0.5, 4.5)) +
      coord_flip(ylim = c(-0.3, 0.3), xlim = c(0.5, 4.5), expand = FALSE, clip = "off") +
      theme(panel.grid = element_blank(),
            panel.background = element_rect(fill = "transparent", colour = "transparent"),
            plot.background = element_rect(fill = "transparent", colour = "transparent"),
            axis.ticks = element_blank(),
            axis.text.x = element_blank(),
            plot.margin = unit(c(0.2, 0.2, 0.2, 0.5),"cm")),
    
    
    # Sickness benefits * Right ideology
    plot_model(lm(Sickness ~ Exp.Group, data = subset(SuppA, Left.Right.Cat == "Right [8-10]"), weights = Weight),
               type = "est",
               title = "",
               colors = "black",
               vline.color = brewer.pal(n = 9, name = "Paired")[5]) +
      geom_hline(yintercept = c(-0.30, -0.25, -0.20, -0.15, -0.10, -0.05, 0.05, 0.10, 0.15, 0.20, 0.25, 0.30), colour = "gray90") +
      geom_point(size = 2.5) +
      geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0, size = 0.7) +
      scale_x_discrete(limits = rev(c("Exp.GroupGroup A:\nGeneral\nproblem\nreminder", 
                                      "Exp.GroupGroup B:\nEconomic\nproblem", 
                                      "Exp.GroupGroup C:\nDeservingness\nproblem", 
                                      "Exp.GroupGroup D:\nEconomic &\ndeservingness\nproblem")),
                       labels = rev(c("Group A: General\nproblem reminder", 
                                      "Group B:\nEconomic problem", 
                                      "Group C:\nDeservingness problem", 
                                      "Group D: Economic &\ndeservingness problem"))) +
      scale_y_continuous(breaks = seq(from = -0.3, to = 0.3, by = 0.05), 
                         labels = c("-0.3", "", "-0.2", "", "-0.1", "", "0.0", "", "0.1", "", "0.2", "", "0.3")) +
      labs(x = NULL, y = NULL, title = NULL) +
      geom_hline(yintercept = c(-0.3, 0.3)) +
      geom_vline(xintercept = c(0.5, 4.5)) +
      coord_flip(ylim = c(-0.3, 0.3), xlim = c(0.5, 4.5), expand = FALSE, clip = "off") +
      theme(panel.grid = element_blank(),
            panel.background = element_rect(fill = "transparent", colour = "transparent"),
            plot.background = element_rect(fill = "transparent", colour = "transparent"),
            axis.ticks.y = element_blank(),
            axis.text.x = element_text(size = 8, vjust = 0.9),
            plot.margin = unit(c(0.2, 0.2, 0.5, 0.5),"cm")),
    
    nrow = 5, heights = c(4.5, 4, 4, 4, 4)),
  
  
  # UNEMPLOYMENT BENEFITS
  ggarrange(
    # Unemployment benefits * Left ideology
    plot_model(lm(Unemployment ~ Exp.Group, data = subset(SuppA, Left.Right.Cat == "Left [0-2]"), weights = Weight),
               type = "est",
               title = "",
               colors = "black",
               vline.color = brewer.pal(n = 9, name = "Paired")[5]) +
      geom_hline(yintercept = c(-0.30, -0.25, -0.20, -0.15, -0.10, -0.05, 0.05, 0.10, 0.15, 0.20, 0.25, 0.30), colour = "gray90") +
      geom_point(size = 2.5) +
      geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0, size = 0.7) +
      scale_y_continuous(breaks = seq(from = -0.3, to = 0.3, by = 0.05), 
                         labels = c("-0.3", "", "-0.2", "", "-0.1", "", "0.0", "", "0.1", "", "0.2", "", "0.3")) +
      labs(x = NULL, y = NULL, title = NULL) +
      geom_rect(ymin = -0.3, ymax = 0.3, xmin = 4.5, xmax = 5.1, colour = "black", fill = "grey90") +
      annotate("text", x = 4.8, y = 0, hjust = 0.5, vjust = 0.5, label = "Unemployment benefits", size = 11/.pt) +
      geom_rect(ymin = 0.3, ymax = 0.4, xmin = 0.5, xmax = 4.5, colour = "black", fill = "grey90") +
      annotate("text", x = 2.5, y = 0.35, hjust = 0.5, vjust = 0.5, label = "Left [0-2]", angle = 270, size = 11/.pt) +
      geom_hline(yintercept = c(-0.3, 0.3)) +
      geom_vline(xintercept = c(0.5, 4.5)) +
      coord_flip(ylim = c(-0.3, 0.4), xlim = c(0.5, 5.0), expand = FALSE, clip = "off") +
      theme(panel.grid = element_blank(),
            panel.background = element_rect(fill = "transparent", colour = "transparent"),
            plot.background = element_rect(fill = "transparent", colour = "transparent"),
            axis.ticks = element_blank(),
            axis.text = element_blank(),
            plot.margin = unit(c(0.5, 0.5, 0.2, 0.2),"cm")),
    
    
    # Unemployment benefits * Left ideology
    plot_model(lm(Unemployment ~ Exp.Group, data = subset(SuppA, Left.Right.Cat == "Centre Left [3-4]"), weights = Weight),
               type = "est",
               title = "",
               colors = "black",
               vline.color = brewer.pal(n = 9, name = "Paired")[5]) +
      geom_hline(yintercept = c(-0.30, -0.25, -0.20, -0.15, -0.10, -0.05, 0.05, 0.10, 0.15, 0.20, 0.25, 0.30), colour = "gray90") +
      geom_point(size = 2.5) +
      geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0, size = 0.7) +
      scale_y_continuous(breaks = seq(from = -0.3, to = 0.3, by = 0.05), 
                         labels = c("-0.3", "", "-0.2", "", "-0.1", "", "0.0", "", "0.1", "", "0.2", "", "0.3")) +
      labs(x = NULL, y = NULL, title = NULL) +
      geom_rect(ymin = 0.3, ymax = 0.4, xmin = 0.5, xmax = 4.5, colour = "black", fill = "grey90") +
      annotate("text", x = 2.5, y = 0.35, hjust = 0.5, vjust = 0.5, label = "Centre left [3-4]", angle = 270, size = 11/.pt) +
      geom_hline(yintercept = c(-0.3, 0.3)) +
      geom_vline(xintercept = c(0.5, 4.5)) +
      coord_flip(ylim = c(-0.3, 0.4), xlim = c(0.5, 4.5), expand = FALSE, clip = "off") +
      theme(panel.grid = element_blank(),
            panel.background = element_rect(fill = "transparent", colour = "transparent"),
            plot.background = element_rect(fill = "transparent", colour = "transparent"),
            axis.ticks = element_blank(),
            axis.text = element_blank(),
            plot.margin = unit(c(0.5, 0.5, 0.2, 0.2),"cm")),
    
    
    # Unemployment benefits * Center ideology
    plot_model(lm(Unemployment ~ Exp.Group, data = subset(SuppA, Left.Right.Cat == "Center [5]"), weights = Weight),
               type = "est",
               title = "",
               colors = "black",
               vline.color = brewer.pal(n = 9, name = "Paired")[5]) +
      geom_hline(yintercept = c(-0.30, -0.25, -0.20, -0.15, -0.10, -0.05, 0.05, 0.10, 0.15, 0.20, 0.25, 0.30), colour = "gray90") +
      geom_point(size = 2.5) +
      geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0, size = 0.7) +
      scale_y_continuous(breaks = seq(from = -0.3, to = 0.3, by = 0.05), 
                         labels = c("-0.3", "", "-0.2", "", "-0.1", "", "0.0", "", "0.1", "", "0.2", "", "0.3")) +
      labs(x = NULL, y = NULL, title = NULL) +
      geom_rect(ymin = 0.3, ymax = 0.4, xmin = 0.5, xmax = 4.5, colour = "black", fill = "grey90") +
      annotate("text", x = 2.5, y = 0.35, hjust = 0.5, vjust = 0.5, label = "Center [5]", angle = 270, size = 11/.pt) +
      geom_hline(yintercept = c(-0.3, 0.3)) +
      geom_vline(xintercept = c(0.5, 4.5)) +
      coord_flip(ylim = c(-0.3, 0.4), xlim = c(0.5, 4.5), expand = FALSE, clip = "off") +
      theme(panel.grid = element_blank(),
            panel.background = element_rect(fill = "transparent", colour = "transparent"),
            plot.background = element_rect(fill = "transparent", colour = "transparent"),
            axis.ticks = element_blank(),
            axis.text = element_blank(),
            plot.margin = unit(c(0.2, 0.5, 0.2, 0.2),"cm")),
    
    
    # Unemployment benefits * Center ideology
    plot_model(lm(Unemployment ~ Exp.Group, data = subset(SuppA, Left.Right.Cat == "Centre Right [6-7]"), weights = Weight),
               type = "est",
               title = "",
               colors = "black",
               vline.color = brewer.pal(n = 9, name = "Paired")[5]) +
      geom_hline(yintercept = c(-0.30, -0.25, -0.20, -0.15, -0.10, -0.05, 0.05, 0.10, 0.15, 0.20, 0.25, 0.30), colour = "gray90") +
      geom_point(size = 2.5) +
      geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0, size = 0.7) +
      scale_y_continuous(breaks = seq(from = -0.3, to = 0.3, by = 0.05), 
                         labels = c("-0.3", "", "-0.2", "", "-0.1", "", "0.0", "", "0.1", "", "0.2", "", "0.3")) +
      labs(x = NULL, y = NULL, title = NULL) +
      geom_rect(ymin = 0.3, ymax = 0.4, xmin = 0.5, xmax = 4.5, colour = "black", fill = "grey90") +
      annotate("text", x = 2.5, y = 0.35, hjust = 0.5, vjust = 0.5, label = "Centre right [6-7]", angle = 270, size = 11/.pt) +
      geom_hline(yintercept = c(-0.3, 0.3)) +
      geom_vline(xintercept = c(0.5, 4.5)) +
      coord_flip(ylim = c(-0.3, 0.4), xlim = c(0.5, 4.5), expand = FALSE, clip = "off") +
      theme(panel.grid = element_blank(),
            panel.background = element_rect(fill = "transparent", colour = "transparent"),
            plot.background = element_rect(fill = "transparent", colour = "transparent"),
            axis.ticks = element_blank(),
            axis.text = element_blank(),
            plot.margin = unit(c(0.2, 0.5, 0.2, 0.2),"cm")),
    
    
    # Unemployment benefits * Right ideology
    plot_model(lm(Unemployment ~ Exp.Group, data = subset(SuppA, Left.Right.Cat == "Right [8-10]"), weights = Weight),
               type = "est",
               title = "",
               colors = "black",
               vline.color = brewer.pal(n = 9, name = "Paired")[5]) +
      geom_hline(yintercept = c(-0.30, -0.25, -0.20, -0.15, -0.10, -0.05, 0.05, 0.10, 0.15, 0.20, 0.25, 0.30), colour = "gray90") +
      geom_point(size = 2.5) +
      geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0, size = 0.7) +
      scale_y_continuous(breaks = seq(from = -0.3, to = 0.3, by = 0.05), 
                         labels = c("-0.3", "", "-0.2", "", "-0.1", "", "0.0", "", "0.1", "", "0.2", "", "0.3")) +
      labs(x = NULL, y = NULL, title = NULL) +
      geom_rect(ymin = 0.3, ymax = 0.4, xmin = 0.5, xmax = 4.5, colour = "black", fill = "grey90") +
      annotate("text", x = 2.5, y = 0.35, hjust = 0.5, vjust = 0.5, label = "Right [8-10]", angle = 270, size = 11/.pt) +
      geom_hline(yintercept = c(-0.3, 0.3)) +
      geom_vline(xintercept = c(0.5, 4.5)) +
      coord_flip(ylim = c(-0.3, 0.4), xlim = c(0.5, 4.5), expand = FALSE, clip = "off") +
      theme(panel.grid = element_blank(),
            panel.background = element_rect(fill = "transparent", colour = "transparent"),
            plot.background = element_rect(fill = "transparent", colour = "transparent"),
            axis.ticks.y = element_blank(),
            axis.text.y = element_blank(),
            axis.text.x = element_text(size = 8, vjust = 0.9),
            plot.margin = unit(c(0.2, 0.5, 0.5, 0.2),"cm")),
    
    nrow = 5, heights = c(4.5, 4, 4, 4, 4)),
  ncol = 2, widths = c(3.16, 2.3)) +
  theme(plot.margin = unit(c(0.1, 0.1, 0.1, 0.1),"cm"),
        panel.background = element_rect(fill = "white", colour = "white"),
        plot.background = element_rect(fill = "white", colour = "white"))



# Saving the output
ggsave(filename = "04-figures and models/FigA1 - Marginal plot - effect vs (categorical) ideology.png", width = 16, height = 22, units = "cm", dpi = 600)






# INTERACTION: AVERAGE EFFECT * INCOME --------------------------------------------------

# Setting 'Group E: Control Group' as a reference category
SuppA$Exp.Group <- factor(SuppA$Exp.Group, 
                          levels = c("Group E:\nControl\nGroup", 
                                     "Group A:\nGeneral\nproblem\nreminder", 
                                     "Group B:\nEconomic\nproblem", 
                                     "Group C:\nDeservingness\nproblem", 
                                     "Group D:\nEconomic &\ndeservingness\nproblem"))


# __ MODELS: Average treatment effects * Income ------------------------------
# Sickness benefits
summary(
  lm(Sickness ~ Exp.Group * Income, 
     data = SuppA, 
     weights = Weight))



# Unemployment benefits
summary(
  lm(Unemployment ~ Exp.Group * Income, 
     data = SuppA, 
     weights = Weight))




# ____ TABLE A6: Regression output -----------------------------------
m5 <- lm(Sickness ~ Exp.Group * Income, data = SuppA, weights = Weight)
m6 <- lm(Unemployment ~ Exp.Group * Income, data = SuppA,  weights = Weight)

stargazer(
  #Interactions with left-right ideology
  m5, m6,
  
  #Specifications
  type = "html", out = "04-figures and models/TabA6 - Income interaction.html",
  intercept.top = TRUE, intercept.bottom = FALSE, no.space = TRUE,
  dep.var.labels = c("Sickness benefits",
                     "Unemployment benefits"),
  covariate.labels = c("Constant",
                       "Group A: General problem reminder",
                       "Group B: Economic problem",
                       "Group C: Deservingness problem",
                       "Group D: Economic and deservingness problem",
                       "Income (group)",
                       "Group A * Income",
                       "Group B * Income",
                       "Group C * Income",
                       "Group D * Income"))


rm(m5, m6)
